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ABSTRACT 


The equation of diffusion of a minor, light gaseous component through an expo- 
nential medium is examined in this detailed report. The equation's one-dimensional 
version may be expressed in dimensionless units in the form 


dn 

dt 


e z 


d 2 n dn 

y z 2 + C 1 +M) JZ +MnJ , 


with /j. the mass ratio of the component particles. With a large variety of initial and 
boundary conditions, the solutions can be expressed in the general form 


n(z, t) 


+ be "' 22 + 


-l+M/2 z 




i=l 


co 

L 


d.(z)e ik - t 


y v £ is the fth zero of the Bessel function J v (x) (v = l and co is the fundamental 
frequency of periodic boundary conditions, if present. 

For t ^ 4 y~\ ~ l,n(z, t) tends, therefore, to approach a steady state regime, 
which may be of periodic character. The former time is also the characteristic time 
scale for transport of the minor component from a source at z = 0 to high altitudes. 
For both hydrogen and helium in the terrestrial atmosphere, this is in the order of 
one day. 

The simple equation is, of course, a most idealized description of diffusion in the 
atmosphere. Nonetheless, the coincidence of the time scales for vertical transport 
to the escape region and the periodic variation found there indicate that a realistic 
evaluation of the density profile of hydrogen or helium should at least, take into ac- 
count the diurnal variation of conditions in the ambient atmosphere. It casts some 
doubt on models of diurnal variation constructed as a superposition of aperiodic sta- 
tionary density profiles. 
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BINARY DIFFUSION IN AN EXPONENTIAL MEDIUM 


by 

Mordehai Liwshitz 
Goddard Space Flight Center 


INTRODUCTION 

Binary diffusion of different gaseous species is common in many physical processes and is a 
typical manifestation of the tendency toward equalization of gross material properties of heteroge- 
neous gaseous mixtures. In this particular case, the properties in question are the relative con- 
centrations of the individual components, their partial pressures, the total pressure of the mixture, 
its temperature, etc. 

Consequently, the gases are set in motion relative to each other by a great variety of general- 
ized forces. Under laboratory conditions of moderate scale with respect to external forces, such 
as gravity, the composition of gaseous mixtures will soon tend toward homogeneity (except near 
the boundaries of containers, where a variety of processes may take place). Diffusion is well 
described by Fick T s second law of diffusion, 


* n io 

<?n 


V ■ ( D i Vn io) • 


( 1 ) 


Here n 10 = rij/N, the ratio of the number density (per cm 3 ) of particles of species 1 to the total 
number density N = (n 1 + n 2 ) /cm 3 for the case of a binary mixture. D x (cm 2 if 1 ) is the diffusion 
coefficient of species 1. In view of the approximate homogeneity soon attained under ordinary con- 
ditions, Equation 1 simplifies to 


dn 10 

~Wt = DV2n • (2) 

This equation, of the same form as the standard equation of heat conduction, is extensively covered 
in References 1 and 2, and admits a great variety of solutions, depending on geometry, boundary 
conditions, and initial conditions. 

Under special conditions in the laboratory and on the grand scale of planetary atmospheres 
and space, Equations 1 and 2 are no longer satisfactory descriptions of the diffusion process. A 
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far more complicated equation, taking into account all of the previously mentioned factors, begins 
with the expression for the velocity of species 1 relative to the other component (Reference 3), 


N 2 


1 10 


n l n 2 ( m 2 m l) Pi P2 , . 

+ Njo lQ gP“ ~£p~ (Ii“£ 2 ) +k T VlogT 


(3) 


Here w a and w 2 (cm/sec) are the velocities of species 1 and 2, nij andm 2 (gm) are their respective 
masses/particle, p(gm/cm t 2 ) is the total specific pressure, p x = n 1 m 1 , p 2 = n 2 m 2 , p = p x + p 2 , 
and Fj and F 2 (cm/sec 2 ) are the accelerations acting on the particles of the two species. K T = D T /D 12 
is the thermal diffusion factor, and T(°K) is the temperature of the gas. 

The time dependent diffusion equation may then be obtained by taking the divergence of the flux 

dn i r -1 

It = " v • [ n i ( ~ x — 2 ) j • (4) 

This equation is not completely rigorous; a more rigorous description is given by the 'Tele- 
grapher’s equation,” including a second derivative with respect to time (Reference 4). For most 
practical purposes, however, the solutions of the two equations are indistinguishable. 

Since the scope of interest lies in the binary diffusion of neutral particles in the atmosphere, 
and gravitational acceleration is independent of mass, the acceleration term in Equation 3 may be 
dropped. The following concerns the diffusion of a minor constituent through an effectively sta- 
tionary major constituent, such that w 2 % 0 , n 2 ~ N . The rigid sphere approximation for the dif- 
fusion coefficient is 


12 


3 

8N o- 
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“ 1 1/2 


aOJ 


yi/2 

“N“ 


(5) 


where k is Boltzmann’s constant, o- is the initial collision cross section, and p = m ly /m 2 . 

To simplify matters further, the one-dimensional diffusion equation will be exclusively con- 
sidered. In view of the above considerations, using Equation 3, and the ideal gas law, (and dropping 
the subscript 1 for n x ), Equation 4 may be written 


dn 

dT 



|“ D (z. t) 


dn (pdN (!-^ +a t) < 9 T( Z , t)\ 

dZ ~ dz ' T(z, t) dz J n 


( 6 ) 


where we assumed that p and a T = K T /n 10 n 20 are approximately constant. Denoting a 2 = - pa 19 and 
a 3 = (l -p + a T ) ap Equation 6 may be obtained in explicit form as 


< 9 n 

Jt 


yi/2 

a i “IT 
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3z 2+ \ N 2 dz + 2NT 1/12 dz) d 2 . + 2N 2 T 1/Z \dz ) \dz ) 
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Changing to X = l/N offers a slight simplification 
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In general, both T and x (and in the worst cases also a x , a 2 , a 3 ) are functions of both z and t ; other 
than numerical solutions of Equation 8 cannot be expected. Though numerical solutions are useful 
and necessary in many situations, they do not provide the desired general insight into the behavior 
of the process investigated. It may, therefore, be profitable to study the properties of an analytic 
solution of a simplified equation. Fortunately the density profile of the main atmosphere is roughly 
exponential, and the temperature, which appears in Equations 7 and 8 in essentially logarithmic 
terms, does not vary over orders of magnitudes. Under the assumption that 


T q , a constant , 
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CO 

P3, 

X Q exp (z/H) , 
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and where 
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bT rt 


Equation 8 simplifies to 
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It is convenient to transform Equation 11 into dimensionless form by changing to dimensionless 
variables 


z 



( 12 ) 


and 


t 1 



(13) 
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Equation 11 becomes 
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-Jt 
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(14) 


where the superscripts on z and t are omitted. It is this equation, and its solution in the region 
from z = 0 to z = oo ? which is discussed in the remaining sections. 

PREVIOUS SOLUTIONS OF THE DIFFUSION EQUATION 
WITH EXPONENTIAL DIFFUSION COEFFICIENT 

Equation 14 has been treated previously by a number of authors (References 5, 6, and 7), 
during investigation of the descent of a pulse of small particles released instantaneously at some 
altitude z 0 in a stationary exponential atmosphere. The solution given in Reference 7 is very use- 
ful for the study of this special problem. But it, and results directly deducible from it, can hardly 
be applied to a realistic description of binary diffusion processes continuing in the upper atmos- 
phere. Such processes are, for instance, the diffusion of hydrogen or helium through an ambient 
medium, composed predominantly of atomic oxygen and molecular nitrogen. 

In changing y = e~ z/2 to a variable, the region (0, °°) in z is transformed to a region (0, 1) 
iny . To avoid complicated boundary condition at y = 1, the authors in Reference 7 now extend the 
region of y to oo (which is equivalent to extending z to -<»), and impose homogeneous boundary condi- 
tions on n for y = 0 and y -<». The initial condition corresponds to an instantaneous particle source 
at time t = 0 and z = z 0 , i.e., n(y, 0)<* 8(y - y 0 ) . 

The resulting solution is 


n(y, t) 


Ht 


y 0 1_M y 1+M exp 



i Im-h 



(15) 


where I, ..is the modified Bessel function of order |/x- l| . Equation 15 predicts the development 
of a constant density profile and its descent with a velocity which is approximately proportional to 
e 2 . Now, this is a reliable description of the particular transient process studied, but fails to 
represent the structure of a component which is supplied by a permanent (though not necessarily 
constant) source located at z Q . Such a source, for instance, supplies neutral atomic hydrogen in 
the vicinity of the turbopause. This becomes most evident for a light component. 

By Equation 3 the net flux of the minor component becomes under the assumptions underlying 
Equation 15, i.e., Equations 9a and 9b 


p _ - e * +fm) 

1 dn (jl 

2y 3y " ^2 n ' 


(16) 
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Choosing for the sake of simplicity the location of the source y 0 at y 0 = 1 (which is permissible 
since the region of y to y-oo has been extended) and writing s = 2yy 0 /t , the instantaneous net flux 
at y = 1 is obtained as 


F(l, s) 


2 exp (-s) I - (s) I 1 - 


s I /3 + 1 ( s ) 
1+ 2 Ig(s) 


(17) 


where /3 = (/x - 1) . The next step is to investigate the asymptotic behavior of f for s -*<»(t-*0), 
and s - 0 (t - oo) . Since the factor multiplying the expression in the fourth set of parentheses in 
Equation 17 is non-negative for all s, it is sufficient to study the asymptotic behavior of the 
expression within the parentheses. 

For t - 0, or s -*oo 


1/3 ( s ) 


( 277 S ) 1/2 


‘ 4/3 2 - 1 

1 " lT8s~ 



(18) 


so that 


s 1/3+1 ( s ) s 3 fj. 
2 I e lTT ~ 2 - 4 + -2 


and the expression within the fourth set of parentheses (Equation 17) tends to 


s s 3 \ _ 1. jx 

V 1 M + 2 2 + 2 _ 4 / " 4 - 2 


Thus for small t the flux will be positive, given n < 1/2 . For t - s - 0 


I /3+l ( s > 


2(/S+ 1) 


(19) 


( 20 ) 


( 21 ) 


and the expression becomes 


[l- M + 0(s)] - 1 - ix . (22) 

For large t the net flux will, therefore, be positive whenever /x < 1, a well known result of 
kinetic theory. The somewhat paradoxical situation is reached that the particle density profile 
(i.e., the bulk of the particles released) moves in one direction, while the net flux is toward the 
opposite direction. Physically this means that the overwhelming majority of particles descends 
with relatively small velocities, while a fast minority soars to high altitudes. 
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But Equation 15 is not suitable for application to the more continuous natural diffusion processes 
for another reason. Since the assumed sources of particles operate over long times, it is presumed 
that some quasi-stationary state is reached. Therefore, a set of boundary conditions and initial 
conditions which provides a solution converging for large t to the solution of the stationary equa- 
tion is a desirable choice (Equation 4). Moreover, the conditions should not display any singu- 
larities, except possibly at the location of the source. 

In order to use the Green’s function, Equation 15, in obtaining the particle distribution ema- 
nating from a continuous source, it is necessary to substitute t - t 0 for t, multiply n(y, t - t 0 ) by 
the specific source yield, s(y 0 . t 0 ), and perform the integration 

n(y, t) = f s(t 0 ) n (y,t- t 0 ) dt 0 . (23) 

Jo 

In view of the complex and singular character of n(y, t - t 0 ), the integral of its product with s (t 0 ) 
may diverge for some values of the parameter 

This divergence, resulting from the superposition of "instantaneous” plane sources, has an 
analog in the theory of heat conduction. 

The superposition of constant point-sources <£(t 0 ) - q yields a temperature 

T(r, t) = 4^7 erfc (^) = (24) 


reducing for t -*<» to 


r ) 47 r «T ’ 


(25) 


where n is the thermal conductivity and erfc (x), is the complimentary error function 


erfc (x) - 1 


/• X 

2 f _ 2 

~ ^ Jo e U 


du . 


T(r ) as given by Equation 25 is well behaved except for the singularity at the origin. 
The corresponding solution for a constant plane source in the plane z = 0 yields 


T(z, t) 


.(a 


1/2 


— z 2 / 4 «t alii 

' r % 


erfc 


2y7Tt 


(26) 


which obviously diverges for large t for all finite z . The difference in behavior between Equa- 
tions 25 and 26 is described to the difference between a singular source of infinitesimal extent and 
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one of infinite extent; this explanation appears to carry over to the case of particle diffusion. At 
any rate, a solution of the type of Equation 26 is not useful since it evidently does not describe the 
physical situation in a realistic manner. 


WELL-BEHAVED SOLUTIONS TO THE BINARY DIFFUSION EQUATION 

For this purpose well-behaved solutions of the diffusion equation solutions are defined as 
non-singular over all t > 0 (except possibly at t = 0), and over all z > 0 (or equivalently 0 < y < 1), 
except possibly at z =0. Moreover, a solution is desirable which for t- 00 converge to the solution 
of the time independent equation 


dn 

dt 


0 


or to a solution of an equation with periodic boundary conditions such as 


(27) 


n(0, t) 


C k e tol 


(28) 


The reason for this choice is an interest in continuous atmospheric diffusion phenomena which 
are likely to display a stationary or periodic behavior. 

Let us then first examine the solution of Equation 27 and Equation 14 under the simple bound- 
ary condition (Equation 28) for large t. In the first case we have to solve the equation of continuity 

( dn \ 

eZ [JI + f jn ) ~ " B ■ (29) 

Here B is the (constant) flux. Its solution is 

n = (1 - M y l j[n 0 (1-/X)-B] e-»“ + Be' 1 } (30) 

where n 0 = n(0). In particular, if B = n 0 (1 -ju) then 


n - n 0 e “ z - (31) 

Solution of Equation 14 with the boundary conditions Equation 28 is very similar to the simplified 
treatment of the earth T s temperature (Reference 9). First we change to a variable y = e“ z/2 , 
which transforms Equation 14 to 


<?y 


dn 

j - (2/J.+ 1) y jp + 4f*n - 


2 

y at 


(32) 
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while the boundary condition becomes 


n(l, t) 


CD 

L 


p ~ik&>t 

H e 


(33) 


Assuming that for t the periodicity of the boundary conditions persists at all altitudes, we look 
for solutions of the form 


n(y, t) = 22 c k u k (y)e ik " t 

n ='-oo 


(34) 


We thus obtain for the kth termu k (y) the equation 


d2u k .. . du k 


y 2 ~ ( 2 M+ l)y + (4^- inoy 2 ) u k 


(35) 


with the boundary condition 


U k (l> = 1 


(36) 


The general solution of Equation 35 is 


U k (y) - y 1+M [a k j| 1 - M | (2i /iKcjy ) + J_ tl _ M| (2i ViKay )] 


(37) 


For a well-behaved solution J_| j_ M) has to be discarded since it diverges at y = 0. Writing? = |k| 
and v - l one obtains therefore 


J ±e ~ a 


±t y 1+M J v ( 2i ) 


(38) 


or, in terms of the Kelvin functions, ber (z) and bei (z) 


u +f (y) - a +? y 1+M |ber v (2 /fe y) + i bei v (2 /fco y)J 

u -£ (y) = a -f y 1+fL e iw ^ber v (lYlco y j - i bei (2-fUo yjj 


(39a) 

(39b) 
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By virtue of the boundary conditions (Equation 36) 


u +£ (y) + u -£ (y) “ 2y 1+/i 1 (2 ber v (2 iJc^j ber v ^2 ita> + bei v (2 too) bei v ^2 ifloo y ) 


(40) 


Here p v (z) = ber v 2 (z) + bei v 2 (z). Writing C ±e = 1/2 (A + iB), the steady state periodic solution of 
Equation 14 is obtained as 


n(z, t ) - C Q e 2 + 


XI ( exp [■ (1 + ^] p ^’ 1 ( 2 *^0 

£= i ' 


ber^ ber v ^2 -flcoe 2/2 j + bei v bei v (2y/£lde z/2 )J 


[A C cos Vcot + Bp sin Pojt ] + 


bei v ^2 iHoXj ber v ^2 ^fTa> e~ z/2 ^ ~ ber v ( 2 l/lcoj bei v ^ 2^fia> e~ 2/2 ^ 
[A g sin loot - Efy cos , 


(41) 


a result valid for o < z < oo . 

The most convenient manner of looking for time dependent solutions of Equation 14 is with the 
use of Laplace transforms 


f( P ) = £ 


i*.>] - r 
J 0 


e~ pt F(t) dt 


(42) 


Changing to y = e 2/2 and applying the Laplace transformation to Equation 14, the equation 

9 d 2 N dN , , x 

y ~j ^2 - (2^+ l)y + 4(/x **■ py 2 ) N - - n(y, + 0)y 2 


(43) 


is obtained for the transform of the density, N(y, p) = £[n(y, t)] . We have not, as yet, 
specified boundary conditions forn(y, t) . For well-behaved solutions it is natural to assume 
thatn- Oasy-O. As for n(l, t) , it will be assumed only that it is some form n x (t) with 
transform Nj(p) . The initial condition n(y, 0+) will not be considered at first because it 
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is seen from the Green's function. Equation 15, that the effect of the initial condition will 
soon die out, even if it is of highly singular character, as long as a finite region of space is 
considered. Moreover, the continuous and more permanent processes usually have a gradual onset, 
represented as a gradual increase of the density from vanishing values at the start. Also, the dis- 
cussion of boundary conditions is greatly facilitated by dealing with the homogeneous part of 
Equation 43. 

The general solution of this equation is 


N(y, P) = y 1+M [aj„ (2ip 1/2 y) +bj. v (2ip 1/2 y)] 


(44) 


(for integral v replace J_„ by Y_„). But for regular solutions in 0 <y < 1, b must vanish. There- 
fore, upon considering the boundary condition N(l, p) = (p), N(y,p)is always of the form 


N(y, p) 


y 1+/i Nj ( P ) 


J„(2ip 1/2 y) 


(45) 


The character, and the mere possibility of evaluation of n(y, t) is therefore crucially dependent 
on the analytic properties of R v - y 1+ ^ J v (2ip 1/2 y |/j v (2ip 1/2 ). Except for the trivial case when 
IM = 1, both numerator and denominator have a branch point at p = 0 in the complex p plane, as 
well as infinite real simple zeros, for v> - 1. In the fraction, however, the branch point is removed, 
and it becomes single valued for all p . It may be written as 


K (P) 


F v (2i P " 2 y) 

F„ (2ip 1/2 ) 


,l + M 


(46) 


where F v (z) are the entire functions of Reference 11. 

When p = pe i0 is written, it is easily seen that the argument 2 p 1/2 e 1 77+0/2 will be real only for 
Q = ±77, i.e,, all (simple) poles of R, (p) will be located on the negative real axis in the p-plane, 
and R^ (p) is a meromorphic complex function with a half -plane of analyticity. This, as will be 
shown in Appendix A, will bear on the existence and properties of the inverse Laplace transform. 

It will, in general, not be possible to invert N(y, p) directly, i.e., to evaluate directly the 
integral 


n(y, t) 


1 
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N(y, p) dp 


(47) 
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and, in most cases, it will be necessary to resort to contour integration. In this context it is also 
necessary to examine the behavior of R v (p) for |p| = p-a o . Writing J u (iz) = (i) v lv(z) , the 
asymptotic expansion of the latter function may be employed to obtain 


iim R, (p) 

p -*oo 


1 im 

p^CD 




0 1/2 / o 1/2 

*2yp y / e 2p y 


]/ 2rrp 1/2 y j ^ 2np 1/2 
= y/x+i/2 i j^e” 2 ^ ( 1 - y ) (cos 6/2 + i sin 0/2 )j 

< y/x+l/2 | j^ e “2jD ( 1 — y ) cos 0 / 2 J , 


< y ^ + V 2 


(48) 


where \o \ < n, y £l. Thus, if N 2 (p) is such that for |p| -a> it tends to zero, another indispensable 
condition for contour integration is met (Appendix). 


It is necessary to examine the behavior of N(y, p) in the half plane of convergence, in order to 
assure that 


N(y, P ) = £[£-»(N)] = N . 


(49) 


Or in other words, in Laplace transformation of the inverse transform n(yj t) = £ -1 [N(y, p)] the 
original transform N(y, p) is recovered, and the resulting inverse transform does not depend on 
the particular abscissa to the right of the abscissa of convergence chosen as the path of integration. 
As outlined in the Appendix, a practically useful sufficient condition is found in Reference 14. 

Assuming N( P )to be analytic in R e (p) >p lr > 0, and assuming it is possible to represent it in 
the form 


N(p) 



(50) 


where (0 <a < 1, e > 0) and where G(p) is bounded in R e (p) >p lr + S > p lr (8 > 0). Then N(p) is the 
Laplace transform of 


r p r +oo 

n( t ) = I e pt N(p) dp . 

Jpr-ico 

As indicated by Equation 45, N(p) = N 1 (p) R v (p) . If then (p) is such that the above condition is 
satisfied, the resulting inverse transform n(t) is proper, regardless of the legitimate manner in 
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which it is evaluated. In particular, since R, (p) is holomorphic in the half plane to the right of the 
first zero of J w (2ip 1/2 ) , N x (p) must be such that the product N(p) fulfills the condition of Equa- 
tion 50. By resorting to contour integration 


n(y, t) 



(51) 


is obtained where is the \th residue of e pt N(y, p). This criterion is applied to a small variety 
of simple, but physically meaningful and instructive, boundary conditions, (t) : 


constant n(l, t) 

n(l. t) = n 0 , 

saturation type 

n(l. t) = n 0 (l-e-«‘) . 

periodic 


uu 

n(l. t) = 


c o ik^t 
e 


M 0 


increasing in t 

n(l, t) = n 0 t a , 

where a > 0 , 


and 


decreasing in t 

n(l, t) = n 0 t 

where 0 < a < 1 . 


Case A. 


N(y, p) 


n 0 R y 

p 


(52a) 


(52b) 


(52c) 


(5 2d) 


(52e) 


(53) 
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The abscissa of convergence is Re (p) = 0. But in view of Equation 50, the representability of 
the inverse n(y, t) as a Laplace transform is not assured, for denoting (y, p) = g( P ) , e = 0 is 
obtained. Since, however, Equation 50 is merely a sufficient condition, a contour integration along 
a suitable path to obtain 


n(y ’ t} = TV 

k=0 


with r k v the kth residue of p* 1 R^ e pt can be performed. Now, N(y, p) has the poles r 0 v at the origin 
and vf (£ / 0) at the zeros of J v (2ip 1/2 ) . These, as was shown, may be written as pj = - y v 2 ? /4 , 
where ±y v ^ are the simple zeros of y v (x). Consequently, 


r 


V 

0 


n 0 y 


2 


since lim R v ( y , p ) 


= y 1+M y v =y 2 - 


Denoting by differentiation with respect to p, also 


(54) 


r e 


v 


n o y 1+M Jv ( 2i P l/2 y) eP ‘ 

pJ v (2ipV2) 


P = 7 v ?e/ 4 


~2n 0 y 1 ^J v (y vgy ) - y ^ t/4 

y v .i Jv+i [y v ,i) e 


(55) 


Thus as the zeros of Bessel functions of consecutive orders interlace, one obtains the alternating 
series 


n(y, t) 





, e Jv+1 {y v ) 


e 






(56) 


which for t -a> converges to the steady state, solution n 0 y 2 , identical with Equation 31. Equation 56 
also obviously satisfies the boundary condition Equation 52a, at all times, i.e., for all t > 0 


limn(y, t) - n Q . (57) 

However, for the convergence at all fixed y of limn(y, t) to vanishing initial value, the t - 0 

t-*o 

theorem based on the conditions of Equation 50 is not applicable since these conditions are not 
satisfied, as indicated above. This criterion may be applied, however, when dealing with boundary 
conditions of type b . 
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Case B. 


N (y. p) 


n 0 „y^j v ( 2ip l/2 y ) 
P(P +«) J v (2 ip 1/2 ) 


(58) 


Equation 58 may be written 


n 0 nR„ (p) 

N(y, P) = • (59) 

Again the abscissa of convergence is the imaginary axis. Then for any S > 0, where S = Re (p) , 
R v/ /(1 +n/p) = G(p) is bounded, and the conditions of Equation 50 are satisfied with & = l > 0. An 
additional term in the series expansion for n(y, t) is obtained, resulting from the pole at p = 

Since « is quite arbitrarily chosen, n / y 2 is assumed, to obviate the evaluation of the residue 
at a second order pole atp = -n is necessary. With this reservation, our series expansion yields 


n(y, t) 



h ( 2" 1/2 y ) e 

hA 1/2 ) 


n t 


00 



«Jv (•/v.ey) 

y v ,i i 4 " rAe) J„ + i (y v ,t) 



(60) 


Again lim (y, t) = n y 2 . But on the strength of the considerations in the Appendix the convergence 

t — * 00 u 

of n(y , t) to zero at t = 0 is assured, since N( y , p) = o(l/p 2 ) for p-°°, and therefore n(y, t) = 0(t) 
for t - 0. For finite t Case A may be approached arbitrarily close by letting 

Case C. 


N(yiP) = R v(yxP) Z](^n^a C k) 


(61) 


again converges for positive real p. The poles are distributed along the negative real axis, and 
along the imaginary axis at p = ±i i spanning all non-negative integers. Using contour integration, 


n(y, t) 


= c o y 2 


. l+p 


CO CD 


{^y v ,i B m -yA K) Jv (y v ,t y) -y v w* 


£ = x 


(y*X + 16m 2 co 2 ) J„ (r v£ ) 


• f + P(y. t) 


(62) 


where P(y, t) is the periodic term in the right member of the steady state solution, Equation 41, 
to which n(y, t) converges for t - <». In view of the convergence criterion Equation 50, convergence 
is also assured for t - 0, if IC k e ika)t is a pure sine series. 
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Case D. 


« 

N (y> p) = n 0 T(a + 1 ) R v (p) p~ (a+l) ( 63 ) 

It is shown in the Appendix that n(y, t ) = o(t a ) for t - <», and, in view of finite solutions for all t , 
this Case can be dismissed. 

Case E, 


N (y> p) - n o r(l-a)R v (p)p 1 a 


(64) 


is shown to vanish as n(y, t) * o(t a ~ 2 ). 


THE EFFECT OF INITIAL CONDITIONS 

Having dealt at length with the dependence of the solution on the boundary conditions, without 
modification of the problem by extension of the investigation to negative values of z (i.e. y > 1), an 
estimate of the effect of initial conditions can be given. The first step is to examine the appropriate 
Green's function G( y , y 0 ; p) obtained from the solution of the equation 

L [ G (y, y 0 ;p)] = s (y~y 0 ) ( 65 ) 

where L is the operator on the left side of Equation 43. Homogeneous boundary conditions at y = 0, 
and y = 1 are chosen. The complete solution n(y, t) will be a sum of the solution to the equation 
with zero initial conditions and the solution with zero boundary conditions. It is easily shown that 


G (y> y 0 ; p) 


ny 

2 v7 


r J v (ay 0 ) Y v (a) , . \ / \ 

" J v (a”) ( a y) + 0 (y o - y) Y v ( a y 0 ) Jv ( a y) + 5 (yyo) Jv ( a y 0 ) Y v ( a y) 


( 66 ) 


with a = 2ip 1/2 . Thus 


N (y. p) = N 0 (y, p) + N. (y, p) . 


( 67 ) 


Here N 0 ( y , P ) denotes the solutions of the previous section, 


N* (y, p) = -4 I" G(y, y 0 ; p) n* (y 0 ) y 0 2 dy„ 

^0 


( 68 ) 
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Two typical distributions may be of interest. 

(A) an initial pulse, as in Reference 7, where n(z, 0 ) = n 0 S(z- z,), or 

n„* y 1 S(y-y 1 ) 

n(y, 0) - — 2 ' 


(69) 


Then from Equations 66 and 68 


N* (y. p) 


nn o yi 1-M y I+/i 


Jv (ay,) Y v (a)J v (ay) 


J„(a) 


- Jv (ay) Y v (ay,) 


where y x < y, and 


™* yj 1 ^ y 1+ ^ 


Y v (a) J v (a yi ) J„ (ay) / x 

■j v (a) -Jv ( a yi) Y V (a-y) 


(70) 


where y > y r The second term on the right side of Equation 69 is identical with the solution in 
Reference 7 while the first term accounts for the adjustment of the initial density to homogeneous 
boundary conditions. 

(B) an initial density corresponding to a previously established steady state n(z, 0) 

= n 0 * e -cri <or> 0 ) , or 

n(y, 0) = n* y 2cr (71) 

Though Equation 71 does not agree with the homogeneous condition stipulated at y = 1, the relation 
of Equation 71 may be assumed to be valid up to y = 1 - e, while beyond that point, it may be assumed 
that n(y , 0) = n 0 * y 2a (y - 1) , and then passes to the limit e - 0. Then 


N (y, p) - 4n 0 Jc( y , yo ; p) y 0 2a+2 dy 0 . (72) 

The integral in Equation 71 is evaluated for two selected values of a: cr l - 1, a 2 = jj. , such that 
N* ( y , p ) may be expressed in terms of a simple combination of Bessel functions. For general a, 
Lommel functions S 2cr v ( z ) and S 2fr _ l i> _ 1 ( z ) would be used. 

(1) a = 1 corresponds to an initial density equal to the steady state density n(z)oc e ~ z . 
Equation 72 yields 


N (y, p) 
4n n 


77y 


l+/i 


Y v ( a ) J v (ay) 


Jv ( a ) 


f y 0 v+ 1 Jv ( a y 0 ) dY o _ Jv ( a y) [ y 0 v+1 Y v ( a y 0 ) dY o ~ Y v ( a y) [ y<T +1 Jv ( a y 0 ) d y 0 

Jo Jy J o J 


7Ty 1+fJ - f r Jv ( a y) 

= Z V-([ Y v( a )Jv + l( a )-Jv( a ) Y v + l(“)] XT^) 


+ 


[Jv ( a y) Y v + x 


(ay) " Y „ (ay) J v+1 



yi+M 3 V (ay) y 2 
Jv ( a ) " a 2 


(73) 
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(2) a = /j. corresponds to diffusive equilibrium, n(t)ae -fiZ , and, in a similar manner, 
Equation 72 yields 


N (y. p) 

4n o 


77V 1+ ^ - , ( a y) y 2 ^ 

25- [Jv (*) Vi <«) - Y v <«) Jv-i <«)] XW " ^T 


yl+M J v (ay) y 2^ 

^2" J„ (a) ” 


( 74 ) 


Substitution y = 0, and y = 1 shows N* (y, P ) to satisfy the boundary conditions, and it is easy to 
verify that Equations 73 and 74 satisfy the inhomogeneous Equation 43. 


THE CHARACTERISTIC STRUCTURE OF SOLUTIONS, 

NUMERICAL RESULTS, AND POSSIBLE APPLICATIONS 

Having obtained a variety of solutions, we now proceed to discuss their salient features. The 
solutions of greater physical interest are Equations 60 and 62. For the former it can be readily 
verified that n(y, 0) = 0. For applying the representation Equation A10 in the Appendix, substi- 
tution z = 2n 1/2 yields 

h (2> t 1/2 y ) - y-> y) 

Jv (2 ^2) " y " " i^y v .t{^-y v 2 x)J v+ i{y^e) ’ (75 


thereby proving the assertion on the convergence of n(y, t) for t - 0. Moreover, by passing to the 
limit n -co 


„ . y 1 Jv(r v ,tyj_ 

y Z-. Jv + i 


y ” 1/2 e ~2(l“y)n 1/2 q 


(76) 


is obtained for y > 1, such that zero initial conditions are also satisfied by Equation 56 for y < 1. 

From our series expansion of n(y, t) it is also evident that the solutions are, indeed, of re- 
laxation type for later times, whatever the initial and boundary conditions. This is in accordance 
with the established aims with respect to the desired solutions. It is also clear that for large t 
(in the order of y~\) the dominant non-periodic, time -dependent term is the first term in the 
series, on account of the factor e 7v (assuming 4^ > yj 2 ^). This term is always of negative 
sign, since J v (y Vtl y )/j v+1 (y Vt < 0 for y < 1. By comparing it with the stationary term, y 2 , a 
characteristic time, t ^ (y) , for approach to steady state conditions can be deduced. For n. » y v l9 
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this becomes 


% (y) 


' V , \ 


1 + log 


' 2 y M-1 3 V {y v .i y ) 

y v .i hn [y v .i) . 


Typical values of are presented in Table 1. 


(77) 


Table 1 

r (y), Characteristic Dimensionless Time for Approach to Steady State Conditions. 


Exponential Depth 


Mass Ratio (^) 


y 

0 

1/2 

1 

0.75 

0.160 

0.180 

0.270 

0.5 

0.385 

0.505 

0.730 

0.25 

0.485 

0.645 

0.935 

0.05 

0.520 

0.685 

1.000 


Table 1 already bears evidence to the easier penetration of lighter particles into a heavier 
medium. The relatively mild dependence of r on p appearing in Table 1 is, however, somewhat 
misleading. By considering dimensional time (using Equations 5 and 13) the dependence is con- 
siderably accentuated. Applying the values in the table to the diffusion of H, 0, N 2 respectively, 
into a medium of o 2 , the respective (y) have to be divided in the ratio V77: Yf: fs . It is as- 
sumed that the cross-section a appearing in Equation 5 does not change appreciably with/x, as well 
as that r 1/32 % r 0 , and r 7/8 * r, . 

On the other hand, the dependence of r ony , appearing in Table 1, explains qualitatively the 
steep slope portion of the density profile progressing in time towards smaller y , i.e., higher al- 
titudes z. r(y ) is, therefore, a rough but quite reliable measure of the diffusion velocity, which in- 
creases noticeably with mounting altitude. 

These results are, of course, not altered substantially if one includes the effect of initial 
conditions. As one can see from Equation 73 and 74 the decay of the initial density is governed by 
the same time dependent series as Equation 55, i.e., 


CD 

E 


J„ (^.gy) 

y v,( Jv+ 1 ,i ) 



As an illustration of results, Figures 1 through 6 represent the evolution of the density profile 
for fj. = 0, 0.5, 1 respectively according to Equations 59 and 60; in the latter, a value of 4^ about 
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Figure 1— Density, n(y, t), versus depth. Figure 2— Density, n(y, t), versus depth, 

y;/x= 0, k = 8.0. y) yL = 1/2, k = 5.5. 
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RELATIVE DENSITY, n(y,t) 



0 0.2 0.4 0.6 0.8 1.0 
EXPONENTIAL DEPTH, y 


Figure 5— Density, n(y, t), versus depth, 

y; m = W, K = °°* 



EXPONENTIAL DEPTH, y 

Figure 6— Density, n(y, t), versus depth, 
y} /m = 1, k - oo. 



EXPONENTIAL DEPTH, y 


Figure 7— Relaxation of density to steady state; 
n(y,t)/n(y, versus y for various times, t. 


halfway between yj* t and y 2 2 has been selected. The 
respective curves show that, even in this case, r as 
computed under the assumption « -cois in the right 
order of magnitude. 

Figure 7 shows the convergence of the density pro- 
file to the steady state for fi = 1/2 . It extends the small 
sample of values presented in Table 1. 

Though it appears that the results of this study 
apply to a very idealized model of the diffusion proc- 
esses going on in the earth T s atmosphere, they con- 
stitute a necessary step in the investigation of more 
complex phenomena. These results were primarily 
intended to serve as zero order solutions of the far 
more complex Equation 8, and provide a convenient 
starting point for a perturbation treatment of this 
equation. 

Nonetheless some meaningful physical information 
may be gleaned by exercising due caution. Of import- 
ance in this respect are the characteristic terms 
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which effectively measure the time of transport for minor components from the source region to 
high altitudes. With the aid of Equations 5 and 13, these are found to be ~10 5 sec. for the diffusion 
of atomic hydrogen starting at ~ 120km and ^3 X 10 s for helium. These light components are sub- 
ject to escape, depending sensitively on the temperature at^500 km, which displays a strong 
diurnal variation. 

In the first approximation the diurnal variation of escape may be regarded as a periodic 
boundary condition, whose period is about the same as the time scale of supply from the source. 

In view of these results (Equation 41), it appears therefore questionable to deduce the density 
profile of minor constituents from the equation of continuity, as has been the usual practice (Ref- 
erences 15, 16 and 17). This will be even more so in the case of hydrogen, where the time of 
depletion, or "life time” with respect to escape (Reference 19) is on the same scale. 
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APPENDIX A 


Some Useful Results from the Theory 
of Laplace Transforms and Complex Function Theory 

Many of the results in this discussion are based on theorems and formulas from the theory 
of Laplace transforms and complex function theory. In fact, the solutions of form F(y,p), are 
obtained directly from solutions of the Laplace transformed equation £[Lf(y, t)] , where 
£(y, t , B/By, B/Bt • • • ) is the differential operator applied to f(y , t ) . 

The immediate question arises whether application of the inverse operator Jl” 1 to F(y, p) yields 
f (y , t) such that £[f(y, t)] = F(y, p) , and for which F(y, p) this can, indeed, be shown to hold. Un- 
fortunately only sufficient conditions are known so far (Reference 14). A rather simple, but readily 
applicable condition is expressed by Theorem 1. 

Let F(p) be analytic in the right half plane Re (p) > p lr > 0, where p = p r + ip. . Let it be possible 
to represent F(p) in this half plane as 


F(P) 


_C_ 

P a 


G(P) 

nl + e 


(Al) 


where (o < a < l, e > o), with G(p) bounded in Re (p) > p lr + S>p lr (S>0) , then F(p) is the Laplace 
transform of the function (t > 0) 


f(t) 


P(277i)" 1 



F(p) dp 


(A2) 


where p r >p lr . Frequent use has been made of this theorem in the examination of "admissible” 
boundary conditions. 

The next question arises with respect to convergence of f(y, t) at t =0. In many cases an 
answer may be found in Theorem 2 (Reference 18). Let F(p) be analytic in Re (p) >p 0r , and let 

F (P> = o (|^r) as p - 00 - 

with /a > 1, independent of the path. Then £“ 1 F(p) converges for t > 0, and represent a function 
f(t) for which 

f(t) = o(t' x ' 1 ) as t - 0 . (A3) 

23 


I 


Similar Abelian theorems (Reference 18) enable the examination of the asymptotic behavior of 
f(t) for t-oo. Theorem 3, below, is a strong theorem of rather wide applicability. Let F(p) be 
analytic in the sector arg(p - p 0 ) <\p (rr/2 < $ < w), except at p 0 . Let F(p) - A(p - p 0 )^ (k real), as 
p-p 0 in this sector. On the rays arg(p-p 0 ) = ±\p let F(p) = o(e n iPt J a s |p| -*». Then 

F(t) ~ AePo 1 as 1 “* 00 * (A4) 

Following the problem of existence of the inverse Laplace transform and its asymptotic be- 
havior for t - 0, and t -<» the question of its evaluation is discussed briefly. Since it is well known 
that this can often be done by simple contour integration, it is sufficient to state two of the neces- 
sary conditions onF(p) (Reference 18), 

(a) F(p) is analytic in the complex p plane except at the poles p 0 , p t * • ■ , all left of K (P)<P 0 r • 

(b) | F( p ) | -Oon the limiting contour C n , i.e., as j p | - 00 . 

It turns out that many of the residues appearing in these results incorporate a factor of the form 

Qv ( w > z ) = 3 P (w)/J v (z) . (A6) 

where w and z are complex variables, connected by the linear relation w = yz. Q v is then a 
meromorphic function of z, and one may invoke the Mittag-Leffler theorem to find a convenient 
representation of (y, z). The branch point of J v (z) at the origin is removed in the quotient, and 
apart from it, J v (z) has simple zeros, yielding the desired representation in the form 

Q* (y- z ) = Qv (y- o) + 22 r k v (f\ + i;) (A7) 

k 

where r ” is the residue of Q v (y, z) at z = z k . For v> -1 all zeros of ] v (z) are real and symmetric 
about the imaginary axis, i.e., z ±k = ±y v k . Also, 

v _ Jv ( y ^>k) 

r+k Jy+l (?v, k ) 

while 

V _ Jy ( y v>k) 

” k Jy+l(^v,k) 


(A8) 


(A9) 


24 



so that 


Q„ (y- z ) 



_1\ Jy ( z kV) 
Z J Jy + l( z k) 


v _ y~< ( 1 _ i 2 \ Jy b'v.k y) 

/ * \ z ^y.k Z ^*^y,k *^y,k / Jy+1 (*yy,k) 


- 2 H 


z2 Jy (■yy.fcy) 


^y.k ( z *-yy%) Jy+l K.k) 


(A10) 


This useful result could have been applied in the inversion of N(p) . This would have entailed, 
however, a complicated examination of the convergence properties of the infinite series in (A10), 
which could be obviated by the use of contour integration. 
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